\defmodule {FDist}

This class provides methods to compute (or approximate)
the distribution functions of special types of
goodness-of-fit test statistics.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bigskip\hrule

\begin{code}
\begin{hide}
/*
 * Class:        FDist
 * Description:  Empirical distributions
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.gof;


public class FDist\begin{hide} {
   private FDist() {}
\end{hide}

   public static double kolmogorovSmirnovPlusJumpOne (int N, double a,
                                                      double x)\begin{hide} {
      final double EPSILONLR = 1.E-15;
      final double EPSILON = 1.0E-290;
      final double NXAPARAM = 6.5;   // frontier: alternating series
      double LogCom;
      double q, p1, q1;
      double Sum = 0.0;
      double term;
      double Njreal;
      double jreal;
      int Sign;
      int j;
      int jmax;

      if (N < 1)
        throw new IllegalArgumentException (
                             "Calling kolmogorovSmirnovPlusJumpOne "+
                             "with N < 1");
      if (a >= 1.0 || a <= 0.0)
        throw new IllegalArgumentException (
                             "Calling kolmogorovSmirnovPlusJumpOne "+
                             "with a outside (0, 1)");
      if (x <= 0.0)
         return 0.0;
      if (x + a >= 1.0)
         return 1.0;
      LogCom = Math.log ((double)N);

      //--------------------------------------------------------------------
      // the alternating series is stable and fast for N*(x + a) very small
      //--------------------------------------------------------------------
      if (N*(x + a) < NXAPARAM && a + x < 0.5) {
         jmax = (int)(N*(x + a));
         for (j = 1; j <= jmax; j++) {
            jreal = j;
            Njreal = N - j;
            q = jreal/N - x;
            if ((q < 0.0 && ((j & 1) != 0)) ||
               ((q > 1.0) && (((N - j - 1) & 1) != 0)))
               Sign = -1;
            else
               Sign = 1;

            // we must avoid log (0.0)
            q1 = Math.abs (q);
            p1 = Math.abs (1.0 - q);
            if (q1 > EPSILON && p1 > EPSILON) {
               term = LogCom + jreal*Math.log (q1) +
                     (Njreal - 1.0)*Math.log (p1);
               Sum += Sign*Math.exp (term);
            }
            LogCom += Math.log (Njreal/(jreal + 1.0));
         }
         // add the term j = 0
         Sum += Math.exp ((N - 1)*Math.log (1.0 + x));
         return Sum*x;
      }

      //---------------------------------------------
      // For N (x + a) >= NxaParam or (a + x) > 0.5,
      // use the non-alternating series.
      //---------------------------------------------

      // EpsilonLR because the distribution has a jump
      jmax = (int)(N*(1.0 - a - x - EPSILONLR));
      for (j = 1; j <= jmax; j++) {
         jreal = j;
         Njreal = N - jreal;
         q = jreal/N + x;
         if (1.0 - q > EPSILON) {
            term = LogCom + (jreal - 1.0)*Math.log (q) + Njreal*Math.log (1.0 - q);
            Sum += Math.exp (term);
         }
         LogCom += Math.log (Njreal/(jreal + 1.0));
      }
      Sum *= x;

      // add the term j = 0
      if (1.0 - x > EPSILON)
         Sum += Math.exp (N*Math.log (1.0 - x));
      return 1.0 - Sum;
   }\end{hide}
\end{code}
 \begin{tabb}
  Similar to
\externalclass{umontreal.iro.lecuyer.probdist}{KolmogorovSmirnovPlusDist}
 but for the case where the distribution function $F$ has a jump of size
 $a$ at a given point $x_0$, is zero at the left of $x_0$,
  and is continuous at the right of $x_0$.
\begin{latexonly}
  The Kolmogorov-Smirnov statistic is defined in that case as
  \eq
    D_N^+(a) = \sup_{a\le u\le 1}
          \left (\hat{F}_N (F^{-1}(u)) - u\right)
             = \max_{\rule{0pt}{7pt} \lfloor 1+aN \le j \le N}
               \left (j/N - F (V_{(j)})\right).
                                    \eqlabel{eq:KSPlusJumpOne}
 \endeq
\iffalse  %%%%%
  and
  \eq
    D_N^-(a) = \sup_{a\le u\le 1} \left (u - \hat F_N (F^{-1}(u))\right)
             = \max_{\rule{0pt}{7pt} \lfloor 1+aN\rfloor \le j \le N}
               \left (F (V_{(j)})-(j-1)/N\right),
  \endeq
 \pierre {It seems that $D_N^-(a)$ has a {\em different\/} distribution
    function. }
\fi  %%%%
  where $V_{(1)},\dots,V_{(N)}$ are the observations sorted by increasing
  order.  The method returns an approximation of
  $P[D_N^+(a) \le x]$ computed via
  \begin{eqnarray}
   P[D_N^+(a) \le x]
    &=& 1 - x \sum_{i=0}^{\lfloor N (1-a-x)\rfloor}
        \binom{N}{i}
        \left (\frac{i}{N} + x \right)^{i-1}
        \left (1 - \frac{i}{N} - x \right)^{N-i}.
        \eqlabel{DistKSJ1} \\[6pt]
    &=& x \sum_{j=0}^{\lfloor N (a+x) \rfloor}
        \binom{N}{j}
        \left (\frac{j}{N} - x \right)^j
        \left (1 - \frac{j}{N} + x \right)^{N-j-1}.
          \eqlabel{DistKSJ2}
  \end{eqnarray}
%  La fonction de r\'epartition est la m\^eme pour  --->  WRONG!
%  \eq
%    D_N^-(a) = \sup_{a\le u\le 1} (u - \hat F (F^{-1}(u)))
%             = \max_{\lfloor 1+aN\rfloor \le j \le N} (F (T_{(j)})-(j-1)/N),
%  \endeq
%  et aussi lorsque le supremum est sur $0\le u \le 1-a$ au lieu de
%  $a\le u\le 1$.
 \hpierre{R\'ef\'erences pour ces formules?}
  The current implementation  uses  formula (\ref{DistKSJ2})
  when $N (x+a) < 6.5$ and $x+a < 0.5$, and uses  (\ref{DistKSJ1})
  when $Nx \ge 6.5$ or $x+a \ge 0.5$.
\end{latexonly}
  Restriction: $0 < a < 1$.
  \end{tabb}
\begin{htmlonly}
   \param{N}{sample size}
   \param{a}{size of the jump}
   \param{x}{positive or negative Kolmogorov-Smirnov statistic}
   \return{the distribution function of the statistic evaluated at \texttt{x}}
\end{htmlonly}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \subsubsection*{Discrete distributions}

\begin{code}

   public static double scan (int N, double d, int m)\begin{hide} {
      return 1.0 - FBar.scan (N, d, m);
   }
}\end{hide}
\end{code}
 \begin{tabb} Returns $F (m)$, the distribution function of the scan statistic
  with parameters $N$ and $d$, evaluated at $m$.
  For a description of this statistic and its distribution, see
  \externalmethod{}{FBar}{scan}{int,double,int},
  which computes its complementary distribution
  $\bar F (m) = 1 - F (m-1)$.
 \end{tabb}
\begin{htmlonly}
   \param{N}{sample size ($\ge 2$)}
   \param{d}{length of the test interval ($\in(0,1)$)}
   \param{m}{scan statistic}
   \return{the distribution function of the statistic evaluated at \texttt{m}}
\end{htmlonly}
